# ==================== 1. 导入依赖库 ====================
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
import tensorflow as tf
from tensorflow.keras.models import Model
from tensorflow.keras.layers import (
    Input, LSTM, Dense, Concatenate, BatchNormalization, Dropout
)
from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau
import warnings
warnings.filterwarnings("ignore")



# 读取表格中的数据
# 1. 读取数据，并初步计算#######################
df1 = pd.read_excel("D:/lxh/北控水务202507/02 事务处理/PAC投加工艺优化/03 数据/data/7月药水1~12.xlsx")
df2 = pd.read_excel("D:/lxh/北控水务202507/02 事务处理/PAC投加工艺优化/03 数据/data/7月药水12~22.xlsx")
df3 = pd.read_excel("D:/lxh/北控水务202507/02 事务处理/PAC投加工艺优化/03 数据/data/7月药水 22~31.xlsx")
df4 = pd.read_excel("D:/lxh/北控水务202507/02 事务处理/PAC投加工艺优化/03 数据/data/2025 08 01.xlsx")
df = pd.concat([df1, df2, df3,df4],axis=0)# 合并

# 提取df中的 时间 一期混合浊 1#沉淀池浊 一期投加量 1#沉淀池水量 PAC配比等列的数据，合并成成新的df
df = df[['时间', '一期混合浊', '一期投加量', '1#沉淀池浊', '1#沉淀池水量', 'PAC配比']]
# 处理变量：把一期混合名字改为进水浊度
df.rename(columns={'一期混合浊': '进水浊度'}, inplace=True)
# 把一期投加量的值除以2
df['一期投加量'] = df['一期投加量'] / 2
# 把'1#沉淀池浊'名字改为出水浊度
df.rename(columns={'1#沉淀池浊': '出水浊度'}, inplace=True)
# 计算PAC单耗：1000*一期投加量*1.23/(1+PAC配比)/1#沉淀池水量
df['PAC单耗'] = 1000 * df['一期投加量'] * 1.23 / (1 + df['PAC配比']) / df['1#沉淀池水量']
# 删除df中的无用列，保留时间、进水浊度、出水浊度、PAC单耗
df = df[['时间', '进水浊度', '出水浊度', 'PAC单耗']]

# 时间 进水浊度 出水浊度 PAC单耗

# 查看数据信息
print(f"数据量：{len(df)} 条，时间范围：{df['时间'].min()} 至 {df['时间'].max()}")
df.head()

# 分解时间列特征
def time_feature_engineering(df):
    """从时间列中提取周期性+长期趋势特征"""
    # 转换为datetime类型（若原始数据不是）
    df["时间"] = pd.to_datetime(df["时间"])
    
    # ------------------------- 周期性特征（解决边界断裂） -------------------------
    # 小时、分钟、星期几
    df["小时"] = df["时间"].dt.hour
    df["分钟"] = df["时间"].dt.minute
    df["星期几"] = df["时间"].dt.weekday  # 0=周一，6=周日
    
    # 正弦/余弦编码（周期24小时、60分钟、7天）
    df["小时_sin"] = np.sin(2 * np.pi * df["小时"] / 24)
    df["小时_cos"] = np.cos(2 * np.pi * df["小时"] / 24)
    df["分钟_sin"] = np.sin(2 * np.pi * df["分钟"] / 60)
    df["分钟_cos"] = np.cos(2 * np.pi * df["分钟"] / 60)
    df["星期几_sin"] = np.sin(2 * np.pi * df["星期几"] / 7)
    df["星期几_cos"] = np.cos(2 * np.pi * df["星期几"] / 7)
    
    # ------------------------- 长期趋势特征（捕捉时间推移） -------------------------
    # 从数据起始时间的总分钟数（反映长期变化）
    start_time = df["时间"].min()
    df["总分钟数"] = (df["时间"] - start_time).dt.total_seconds() // 60  # 整数分钟
    
    # ------------------------- 特征筛选 -------------------------
    time_features = [
        "小时_sin", "小时_cos", "分钟_sin", "分钟_cos", 
        "星期几_sin", "星期几_cos", "总分钟数"
    ]
    core_features = ["进水浊度", "出水浊度", "PAC单耗"]  # 核心业务特征
    
    return df, time_features, core_features

# 执行特征工程
df, time_features, core_features = time_feature_engineering(df)

def data_cleaning_and_scaling(df, time_features, core_features):
    """数据清洗与标准化"""
    # ------------------------- 缺失值处理 -------------------------
    # 检查缺失值（模拟数据无缺失，真实数据需处理）
    missing = df.isnull().sum()
    if missing.sum() > 0:
        print(f"缺失值统计：{missing[missing>0]}")
        # 线性插值补全（适用于短时间缺失）
        df = df.interpolate(method="linear", limit_direction="both")
    
    # ------------------------- 异常值处理（基于IQR） -------------------------
    def remove_outliers(series, threshold=3):
        """基于IQR剔除异常值"""
        q1 = series.quantile(0.25)
        q3 = series.quantile(0.75)
        iqr = q3 - q1
        lower = q1 - threshold*iqr
        upper = q3 + threshold*iqr
        return series.clip(lower, upper)
    
    # 对核心特征（浊度、投加量）去异常
    for col in core_features:
        df[col] = remove_outliers(df[col])
    
    # ------------------------- 特征标准化 -------------------------
    # 时间特征标准化（范围[-1,1]，与正弦/余弦编码兼容）
    scaler_time = MinMaxScaler(feature_range=(-1, 1))
    df[time_features] = scaler_time.fit_transform(df[time_features])
    
    # 核心特征标准化（范围[0,1]，因浊度/投加量为非负）
    scaler_core = MinMaxScaler(feature_range=(0, 1))
    df[core_features] = scaler_core.fit_transform(df[core_features])
    
    return df, scaler_time, scaler_core

# 执行清洗与标准化
df, scaler_time, scaler_core = data_cleaning_and_scaling(df, time_features, core_features)
df.head()




def create_sequences(df, time_features, core_features, time_steps=144, pred_steps=1):
    """
    构造LSTM输入序列和目标值
    :param time_steps: 历史窗口长度（10分钟×time_steps）
    :param pred_steps: 预测未来步数（1=10分钟，2=20分钟...）
    """
    X_time, X_core, y = [], [], []  # 时间特征序列、核心特征序列、目标值
    
    # 总样本数 = 数据量 - 窗口长度 - 预测步数 + 1
    total_samples = len(df) - time_steps - pred_steps + 1
    print(f"构造样本数：{total_samples}")
    
    for i in range(total_samples):
        # 历史窗口（时间特征+核心特征）
        start_idx = i
        end_idx = i + time_steps
        
        # 时间特征序列（形状：(time_steps, len(time_features))）
        seq_time = df[time_features].iloc[start_idx:end_idx].values
        X_time.append(seq_time)
        
        # 核心特征序列（形状：(time_steps, len(core_features))）
        seq_core = df[core_features].iloc[start_idx:end_idx].values
        X_core.append(seq_core)
        
        # 目标值（未来pred_steps步的出水浊度和投加量）
        # 这里预测未来1步的出水浊度（第core_features[1]列）和投加量（第core_features[2]列）
        target = df[[core_features[1], core_features[2]]].iloc[end_idx:end_idx+pred_steps].values
        y.append(target.flatten())  # 展平为一维数组（2个值）
    
    # 转换为numpy数组
    X_time = np.array(X_time)  # 形状：(样本数, time_steps, len(time_features))
    X_core = np.array(X_core)  # 形状：(样本数, time_steps, len(core_features))
    y = np.array(y)            # 形状：(样本数, 2)
    
    return X_time, X_core, y

# 参数设置（时间步长=144=1天，预测未来1步=10分钟）
TIME_STEPS = 144
PRED_STEPS = 1

# 构造样本
X_time, X_core, y = create_sequences(df, time_features, core_features, TIME_STEPS, PRED_STEPS)
print(f"X_time形状：{X_time.shape}, X_core形状：{X_core.shape}, y形状：{y.shape}")

def time_series_split(X_time, X_core, y, train_ratio=0.7, val_ratio=0.2, test_ratio=0.1):
    """按时间顺序划分数据集"""
    total_samples = len(X_time)
    train_end = int(total_samples * train_ratio)
    val_end = train_end + int(total_samples * val_ratio)
    
    # 训练集：前70%
    X_time_train, X_core_train, y_train = X_time[:train_end], X_core[:train_end], y[:train_end]
    # 验证集：中间20%
    X_time_val, X_core_val, y_val = X_time[train_end:val_end], X_core[train_end:val_end], y[train_end:val_end]
    # 测试集：最后10%
    X_time_test, X_core_test, y_test = X_time[val_end:], X_core[val_end:], y[val_end:]
    
    return (X_time_train, X_core_train, y_train), (X_time_val, X_core_val, y_val), (X_time_test, X_core_test, y_test)

# 执行划分
(train_time, train_core, train_y), (val_time, val_core, val_y), (test_time, test_core, test_y) = time_series_split(
    X_time, X_core, y, train_ratio=0.7, val_ratio=0.2, test_ratio=0.1)

print(f"训练集样本数：{len(train_time)}")
print(f"验证集样本数：{len(val_time)}")
print(f"测试集样本数：{len(test_time)}")

def build_lstm_model(time_steps, time_feats_dim, core_feats_dim):
    """构建多输入LSTM模型"""
    # ------------------------- 时间特征输入分支 -------------------------
    time_input = Input(shape=(time_steps, time_feats_dim), name="time_input")
    time_lstm = LSTM(64, return_sequences=False, dropout=0.2)(time_input)  # 提取时间模式
    time_dense = Dense(32, activation="relu")(time_lstm)  # 压缩时间特征
    
    # ------------------------- 核心特征输入分支 -------------------------
    core_input = Input(shape=(time_steps, core_feats_dim), name="core_input")
    core_lstm = LSTM(128, return_sequences=False, dropout=0.2)(core_input)  # 提取业务模式
    core_dense = Dense(64, activation="relu")(core_lstm)  # 压缩核心特征
    
    # ------------------------- 特征融合与输出 -------------------------
    merged = Concatenate()([time_dense, core_dense])  # 融合时间与业务特征
    outputs = Dense(2, activation="linear", name="output")(merged)  # 预测出水浊度和投加量
    
    # 定义模型
    model = Model(inputs=[time_input, core_input], outputs=outputs)
    model.compile(
        optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),
        loss="mse",  # 均方误差（回归任务）
        metrics=["mae"]  # 平均绝对误差
    )
    
    return model

# 模型参数
TIME_FEATS_DIM = len(time_features)  # 时间特征维度（7）
CORE_FEATS_DIM = len(core_features)  # 核心特征维度（3）

# 构建模型
model = build_lstm_model(TIME_STEPS, TIME_FEATS_DIM, CORE_FEATS_DIM)
model.summary()  # 打印模型结构

# 回调函数（早停+学习率衰减）
callbacks = [
    EarlyStopping(monitor="val_loss", patience=15, restore_best_weights=True),  # 早停
    ReduceLROnPlateau(monitor="val_loss", factor=0.5, patience=5, min_lr=1e-6)  # 学习率衰减
]

# 训练模型
history = model.fit(
    x={"time_input": train_time, "core_input": train_core},  # 多输入
    y=train_y,
    validation_data=(
        {"time_input": val_time, "core_input": val_core}, 
        val_y
    ),
    batch_size=64,  # 批次大小（根据GPU内存调整）
    epochs=100,     # 最大轮次（早停会提前终止）
    callbacks=callbacks,
    verbose=1
)

# 绘制训练曲线
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(history.history["loss"], label="训练损失")
plt.plot(history.history["val_loss"], label="验证损失")
plt.title("损失曲线")
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(history.history["mae"], label="训练MAE")
plt.plot(history.history["val_mae"], label="验证MAE")
plt.title("MAE曲线")
plt.legend()
plt.tight_layout()
plt.show()

# 模型评估（测试集）


# 预测测试集
test_pred = model.predict({"time_input": test_time, "core_input": test_core})
test_pred = scaler_core.inverse_transform(test_pred)  # 反归一化（核心特征）
test_true = scaler_core.inverse_transform(test_y)     # 反归一化（真实值）

# 计算指标
mse = mean_squared_error(test_true, test_pred)
mae = mean_absolute_error(test_true, test_pred)
r2 = r2_score(test_true, test_pred)

print(f"测试集评估结果：")
print(f"MSE（出水浊度+投加量）: {mse:.2f}")
print(f"MAE（出水浊度+投加量）: {mae:.2f}")
print(f"R²分数: {r2:.2f}")

# 可视化预测对比（以出水浊度为例）
plt.figure(figsize=(12, 6))
plt.plot(test_true[:, 0], label="真实出水浊度", alpha=0.6)
plt.plot(test_pred[:, 0], label="预测出水浊度", alpha=0.6)
plt.title("测试集出水浊度预测对比")
plt.xlabel("样本索引")
plt.ylabel("出水浊度(NTU)")
plt.legend()
plt.show()



# 预测0-3200NTU进水时的PAC单耗和出水浊度
def predict_target_scenario(model, df, time_features, core_features, scaler_core, 
                            target_inlet_range=(0, 3200), time_steps=144, pred_steps=1):
    """
    预测进水浊度在指定范围内的出水浊度和投加量
    :param target_inlet_range: 目标进水浊度范围（NTU）
    """
    # ------------------------- 筛选符合条件的历史窗口 -------------------------
    # 核心特征中的进水浊度列索引（core_features[0]）
    inlet_col_idx = core_features.index("进水浊度(NTU)")
    
    # 遍历数据，找到进水浊度在目标范围内的窗口
    valid_samples = []
    for i in range(len(df) - time_steps - pred_steps + 1):
        window_inlet = df[core_features[inlet_col_idx]].iloc[i:i+time_steps]
        if (window_inlet.min() >= target_inlet_range[0]) and (window_inlet.max() <= target_inlet_range[1]):
            # 提取时间特征序列和核心特征序列
            seq_time = df[time_features].iloc[i:i+time_steps].values.reshape(1, time_steps, -1)
            seq_core = df[core_features].iloc[i:i+time_steps].values.reshape(1, time_steps, -1)
            valid_samples.append((seq_time, seq_core))
    
    # 转换为数组
    valid_time = np.concatenate([s[0] for s in valid_samples])
    valid_core = np.concatenate([s[1] for s in valid_samples])
    
    # ------------------------- 预测 -------------------------
    pred = model.predict({"time_input": valid_time, "core_input": valid_core})
    pred = scaler_core.inverse_transform(pred)  # 反归一化
    
    # 提取进水浊度（取窗口内的平均值作为代表）
    window_inlet_mean = [df[core_features[inlet_col_idx]].iloc[i:i+time_steps].mean() for i in range(len(valid_samples))]
    
    return window_inlet_mean, pred

# 执行预测（进水浊度0-3200NTU）
target_inlet_min, target_inlet_max = 0, 3200
window_inlet_mean, pred = predict_target_scenario(
    model, df, time_features, core_features, scaler_core,
    target_inlet_range=(target_inlet_min, target_inlet_max),
    time_steps=TIME_STEPS
)

# 提取出水浊度和投加量（pred的第二列是投加量）
pred_outlet = pred[:, 0]
pred_dosing = pred[:, 1]

# 绘制预测曲线
# 可视化进水浊度与出水浊度/投加量的关系
plt.figure(figsize=(14, 6))

# 出水浊度 vs 进水浊度
plt.subplot(1, 2, 1)
plt.scatter(window_inlet_mean, pred_outlet, alpha=0.6, color="blue")
plt.plot([target_inlet_min, target_inlet_max], [target_inlet_min, target_inlet_max], 
         "r--", label="理想情况（出水=进水）")  # 理论下限（无处理）
plt.title("进水浊度 vs 预测出水浊度")
plt.xlabel("进水浊度(NTU)")
plt.ylabel("预测出水浊度(NTU)")
plt.legend()
plt.grid(True)

# 投加量 vs 进水浊度
plt.subplot(1, 2, 2)
plt.scatter(window_inlet_mean, pred_dosing, alpha=0.6, color="red")
plt.title("进水浊度 vs 预测投加量")
plt.xlabel("进水浊度(NTU)")
plt.ylabel("预测投加量(kg/h)")
plt.grid(True)

plt.tight_layout()
plt.show()